#include <iostream>
#include <vector>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

// Function to perform discrete least square via normal equations
vector<double> discreteLeastSquare(const vector<double>& x, const vector<double>& y) {
    int N = x.size();
    MatrixXd A(N, 3);
    VectorXd b(N);

    // Constructing matrix A and vector b
    for (int i = 0; i < N; ++i) {
        A(i, 0) = 1;
        A(i, 1) = x[i];
        A(i, 2) = x[i] * x[i];
        b(i) = y[i];
    }

    // Solving the normal equations
    VectorXd coeffs = (A.transpose() * A).ldlt().solve(A.transpose() * b);
    return {coeffs[0], coeffs[1], coeffs[2]};
}

int main() {
    // Input data
    vector<double> x = {0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0};
    vector<double> y = {2.9, 2.7, 4.8, 5.3, 7.1, 7.6, 7.7, 7.6, 9.4, 9.0, 9.6, 10.0, 10.2, 9.7, 8.3, 8.4, 9.0, 8.3, 6.6, 6.7, 4.1};

    // Compute the coefficients
    vector<double> coeffs = discreteLeastSquare(x, y);

    // Output the coefficients
    cout << "The coefficients are: " << endl;
    cout << "a0: " << coeffs[0] << endl;
    cout << "a1: " << coeffs[1] << endl;
    cout << "a2: " << coeffs[2] << endl;

    return 0;
}
